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Mixed-Mode Decohesion Finite Elements for the 
Simulation of Delamination in Composite 

Materials 

P.P. Camanho* and C.G. Davila^ 


Abstract 

A new decohesion element with mixed-mode capability is proposed and 
demonstrated. The element is used at the interface between solid finite ele- 
ments to model the initiation and non-self-similar growth of delaminations. 
A single relative displacement-based damage parameter is applied in a soft- 
ening law to track the damage state of the interface and to prevent the 
restoration of the cohesive state during unloading. The softening law for 
mixed- mode delamination propagation can be applied to any mode interac- 
tion criterion such as the two-parameter power law or the three-parameter 
Benzeggagh-Kenane criterion. To demonstrate the accuracy of the predic- 
tions and the irreversibility capability of the constitutive law, steady-state 
delamination growth is simulated for quasi-static loading-unloading cycles 
of various single mode and mixed-mode delamination test specimens. 


1 Introduction 

Interlaminar damage (delamination) is one of the predominant forms of failure in 
many laminated composites systems, especially when there is no reinforcement in 
the thickness direction. Delamination as a result of impact or a manufacturing 
defect can cause a significant reduction in the compressive load-carrying capacity 
of a structure. The stress gradients that occur near geometric discontinuities such 
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as ply drop-offs, stiffener terminations and flanges, bonded and bolted joints, and 
access holes promote delamination initiation, trigger intraply damage mechanisms, 
and may cause a significant loss of structural integrity. 

The fracture process in high performance composite laminates is quite complex, 
involving not only delamination, but also intralaminar damage mechanisms (e.g. 
transverse matrix cracking, fiber fracture). For effective predictive capabilities, 
failure analysis tools for the different failure modes are required. 

The simulation of delamination in composites is usually divided into delami- 
nation initiation and delamination propagation. Delamination initiation analyses 
are usually based on stresses and use criteria such as the quadratic interaction of 
the interlaminar stresses in conjunction with a characteristic distance [1],[2]. The 
characteristic distance is an averaging length that is a function of geometry and 
material properties, so its determination always require extensive testing. 

Most analyses of delamination growth apply a fracture mechanics approach, 
and evaluate energy release rates G for self-similar delamination growth [3]- [9]. 
The energy release rates are usually evaluated using the virtual crack closure tech- 
nique (VCCT) proposed by Rybicki and Kanninen [10]. The VCCT technique is 
based on Irwin’s assumption that when a crack extends by a small amount, the 
energy released in the process is equal to the work required to close the crack to 
its original length. The Mode I, Mode II, and Mode ITT energy release rates, Gj, 
Gn and Gni respectively, can then be computed from the nodal forces and dis- 
placements obtained from the solution of a finite element model. The approach is 
computationally effective since the energy release rates can be obtained from only 
one analysis. Although valuable information concerning the onset and the stabil- 
ity of delamination can be obtained using the VCCT, its use in the simulation of 
delamination growth may require complex moving mesh techniques to advance the 
crack front when the local energy release rates reach a critical value [11]. Further- 
more, an initial delamination must be defined and, for certain geometries and load 
cases, the location of the delamination front might be difficult to determine. 

The use of decohesion elements placed at the interfaces between laminae can 
overcome some of the above difficulties. Decohesion elements are based on a 
Dudgale-Barenblatt cohesive zone approach [12], [13], which can be related to 
Griffith’s theory of fracture when the cohesive zone size is negligible when com- 
pared with characteristic dimensions, regardless of the shape of the constitutive 
equation [14]. These elements use failure criteria that combine aspects of strength- 
based analysis to predict the onset of the softening process at the interface and 
Fracture Mechanics to predict delamination propagation. A main advantage of the 
use of decohesion elements is the capability to predict both onset and propagation 
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of delamination without previous knowledge of the crack location and propagation 
direction. 

Decohesion elements can be divided into two main groups: continuous interface 
elements and point decohesion elements. Different types of continuous decohesion 
elements have been proposed, ranging from zero-thickness volumetric elements con- 
necting solid elements [15], finite-thickness volumetric elements connecting shell 
elements [16], and line elements [17]- [18]. Point decohesion elements are identical 
to non-linear spring elements connecting nodes [19], [20]. A common feature of 
previously developed decohesion elements is the absence of an interaction crite- 
rion for the prediction of softening onset under mixed-mode loading and the use 
of simplified interaction criteria of the energy release rates for the prediction of 
delamination propagation. However, experimental evidence shows that for some 
resins (e.g. epoxies) the dependence of the fracture toughness on the mode ratio 
may not be expressed by a simplified expression [21]. Under mixed- mode loading 
conditions, the propagation of delamination should be predicted using physically 
sound criteria. 

The objective of the current work is to develop zero-thickness volumetric deco- 
hesion elements able to capture delamination onset and growth under mixed-mode 
loading conditions. A quadratic interaction between the tractions is proposed 
to predict softening onset. A criterion able to capture the mixed-mode fracture 
toughness under different mode ratios is used to predict delamination propagation. 
The capabilities of the proposed formulation are assessed simulating double can- 
tilever beam (DCB), end- notch flexure (ENF) and mixed- mode bending (MMB) 
test specimens, and comparing the predictions with experimental data. 


2 Decohesion element formulation 

2.1 Element kinematics 

The zero- thickness decohesion elements with 8-nodes (shown in Figure 1) and 18- 
nodes are proposed to simulate the resin-rich layer connecting the several laminae 
of a composite laminate. The constitutive equation of zero- thickness decohesion 
elements is established in terms of relative displacements (also called displacement 
discontinuities [22]) and tractions across the interface. 

The definition of the relative displacements for an element with a general ori- 
entation in space is obtained using a procedure based on the work of Ahmad [23] 
and Beer [24]. The vector defining the relative displacement in global coordinates, 
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Figure 1: Zero-thickness decohesion element 
A, can be obtained as: 


A* NkUfa NkUfci (1) 

where u ^ are the displacements in the i direction of the k top and bottom 
nodes of the element, respectively. are standard Lagrangian shape functions. 

For a general element shape and alignment, the normal and tangential relative 
displacements must be determined in local coordinates. The tangential plane at a 
given point is spanned by two vectors, v^ and v^, obtained by differentiating the 
global position vector with respect to the natural (local) coordinates: 

% = X i,Z’ V % = X i,V ( 2 ) 

Defining an isoparametric element, the global position vector is obtained as: 


From (2) and (3): 


^k^ki 


( 3 ) 


{J^k^ki) ,£ Nk£%ki ( 4 ) 

' % 

{Nk%ki) yrj ^k,r]^ki (^) 

Although and are, generally, not orthogonal to each other, their vector 
product defines a surface normal. Therefore, the local normal coordinate vector is 
obtained as: 


Vn = (V£ X V^) ||V£ X vj (6) 

The tangential coordinates are then obtained as: 

V s = v e ||v^|| _1 
Vt = v„ X v s 


4 


( 7 ) 

(8) 



The components of v n , v s , and v t represent the direction cosines of the local 
coordinate system to the global coordinate system, thus defining the transforma- 
tion tensor @ si . Using (1), the relative displacements can then be obtained in local 
coordinates as: 


®si^i ® si^ k^ki ^sik^ki (9) 

The constitutive operator of the decohesion element, £ sr , relates the element 
tractions, r s , to the element relative displacements, 6 r : 

r s =S sr D sr 8 r (10) 

where 6 sr is the Kronecker delta. 

The coefficients D sr of the element constitutive operator can be used to simulate 
elaborate mechanical behaviors, including the mechanics of interfacial decohesion 
and crack propagation, and will be discussed later 1 . An important characteris- 
tic of the proposed method is that, unlike thin continuum elements (degenerate 
continuum elements), the stiffness of the interface before softening onset, referred 
here as the penalty stiffness , is not a function of the discretization, but is defined 
by the coefficients D sr . Some authors [25] have proposed the definition of the 
penalty stiffness as a function of the interface thickness, £, and elastic moduli of 
the interface (£ 3 , G 13 and G 23 ) as: £33 = £ 3 /^, £11 = 2G\$/t, £22 = ‘ZGm/t. 

The decohesion element stiffness matrix and internal load vector can be ob- 
tained from the principle of virtual work: 


j d8 s r s dT - fkidu ki = 0 

From (9), and considering a geometrically linear problem: 


(ii) 


£si&Tsdr fki 0 


( 12 ) 


The first term of (12) represents the decohesion element internal load vector. 
From (9) and (10): 


rvzd J T''Uzv fki 

(13) 

K-ikvz'U'zv fki 

(14) 


^or simple contact elements D sr are the penalty parameters: D sr — 0 if 63 > 0 , and D33 = K 
if <S 3 < 0. 
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The decohesion element stiffness matrix is then: 


H-l /*+ 1 


Kikzv / B sik 6 S r D sr B rvz dT 

Jr 


f - 1 J-i 


B sik 6 sr D sr B rvz || v € x v,j|| (15) 


Care must be exercised in the choice of the integration scheme used to ob- 
tain the stiffness matrix and the nodal force vector. Several investigations have 
shown the superiority of using Newton-Cotes integration techniques over tra- 
ditional Gaussian integration techniques in decohesion elements [26]- [30]. Us- 
ing eigenmode analysis of the element stiffness matrices it has been shown that 
Gaussian integration can cause undesired spurious oscillations of the traction field 
when large traction gradients are present over an element [28] , [29] . 

Another relevant issue related with the integration scheme is the number of 
integration points used. Analyses of problems involving crack propagation and 
softening behavior have shown that the use of full integration was superior to 
the use of reduced integration schemes [31]. However, Alfano and Crisfield [32] 
have shown that for linear 4-node decohesion elements increasing the number of 
integration points from 2 (corresponding to full integration) to 20 results in an 
increase of spurious oscillations in the load-displacement curve and, consequently, 
in a less robust solution algorithm. For the above reasons, Newton-Cotes full 
integration is used in the decohesion element proposed here. 


2.2 Proposed constitutive equation 

2.2.1 Single-mode delamination 

The need for an appropriate constitutive equation in the formulation of the decohe- 
sion element is fundamental for an accurate simulation of the interlaminar cracking 
process. It is considered that there is a process zone or cohesive zone ahead of 
the delamination tip. Figure 2 represents the cohesive zone in specimens loaded in 
pure Mode II (Figure 2-a)) and in pure Mode I (Figure 2-b)). Figure 2 also illus- 
trates the constitutive behaviour for pure Mode I, pure Mode II, and pure Mode 
III loading. The concept of cohesive zone was initially proposed by Barenblatt [13] 
and using such a concept the singularity at the crack tip is removed. 

Physically, the cohesive zone represents the coalescence of crazes in the resin 
rich layer located at the delamination tip and reflects the way by which the material 
loses load-carrying capacity [33]. Ungsuwarungsri and Knauss [33] considered that 
if the size of the process zone is narrow compared to the size of the specimen, a 
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softening material behavior confined to a thin layer adjacent to the crack plane is 
a realistic scheme for the simulation of crack growth. Needleman [34] considered 
that cohesive zone models are particularly attractive when interfacial strengths 
are relatively weak when compared with the adjoining material, as is the case in 
composite laminates. 




a) Mode II or Mode III 


b) Mode I 


Figure 2: Pure mode constitutive equations 

It has been shown that cohesive zone approaches can be related to Griffith’s 
theory of fracture if the area under the traction-relative displacement relation is 
equal to the corresponding fracture toughness [35] (see Figure 2), regardless of its 
shape. Furthermore, Crisfield [32] has shown that when the relative displacements 
and Si shown in Figure 2 are coincident (corresponding to a sudden load drop 
to zero) a perfectly brittle fracture is simulated. A model for brittle fracture must 
be able to capture the high stress gradients at the crack tip with sufficiently fine 
mesh densities or singular elements. 

For pure Mode I and pure Mode II or Mode III loading the bi- linear soften- 
ing constitutive behaviour represented in Figure 2 is used. A high initial stiffness 
(penalty stiffness , K) is used to hold the top and bottom faces of the decohesion 
element together in the linear elastic range (point 1 in Figure 2). For pure Mode I, 
II or III loading, after the interfacial normal or shear tractions attain their respec- 
tive interlaminar tensile or shear strengths (point 2 in Figure 2), the stiffnesses are 
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gradually reduced to zero. The onset displacements are obtained as: £3 = N/K , 
62 = S/K and 6° = T/K, where N is the interlaminar tensile strength, and S and 
T are the interlaminar shear strengths. 

The area under the traction-relative displacement curves is the respective (Mode 
I, II or III) fracture toughness (G/c, Gnc and Gmc respectively) and defines the 
final relative displacements, 63, and <5*(, corresponding to complete decohesion: 


G IC (16) 

Gnc (17) 

Gmc (18) 

The final displacements are then obtained as: = 2Gjc/N , 8{ = 2 Gjjc/S 

and S{ = 2Gmc/T. 

Once a crack is unable to transfer any further load (point 5 in Figure 2), all 
the penalty stiffnesses revert to zero. However, it is necessary to avoid the inter- 
penetration of the crack faces. The contact problem is addressed by re-applying 
the normal penalty stiffness when interpenetration is detected. 

In order to formulate the complete constitutive equation, the unloading behav- 
ior must be defined. It is considered that a softening point unloads towards the 
origin, as shown in Figure 2. Using the following operator: 

( 0 <£= x < 0 

<*> = \ (19) 

y x <= x > 0 

the loading condition can be formulated in terms of a state variable defined as the 
maximum relative displacement , 6 max , suffered by the point: 

Mode II or III : <5™* = max {<5™“, |^|} , * = 1,2 (20) 

Mode I : Sg** = max {<$”**, M , with 6” ax >0 (21) 



and using a loading function, F , defined as: 



Mode II or III : 

f( \6i\-er*] 

II 

cmax\ 

* 2 = 12 
cmax 1 b L 1 * 

( 22 ) 

Mode I : 

F(<§ 3 - $T) 

1 1 

CC CO 

II 

8 g““) 

cmax J 

®3 

(23) 


with 6^ > 0. 


Using <5 max in the constitutive equation, the irreversibility of damage is taken 
into account. This is shown in Figure 2: if the relative displacement decreases, the 
point unloads elastically towards the origin with a reduced, secant, stiffness (point 
3 in Figure 2). 

The irreversible, bi-linear, softening constitutive behaviour shown in Figure 2 
have been developed in previous work [17], [32], [36], and can be defined as: 


Ti 


di 


r K8i <*= 6™ x <8° 

< (1 -di)K6i <= 8°< <5f ax < S{ 
0 4= 6™°* > S{ 


8j_ Z 6 J1 

eriti-fii)* 


1,2,3; di e [0,1] 


(24) 


(25) 


In order to avoid interpenetration of the crack faces, the following condition is 
introduced: 


r 3 = KSs <= 63 < 0 (26) 

The properties required to define the interfacial behavior are the penalty stiff- 
ness, K, the corresponding fracture toughness, Gnc and Gnic , and the 

corresponding interlaminar normal tensile or shear strengths, iV, S or T respec- 
tively. 


2.2.2 Mixed-mode delamination 

In structural applications of composites, delamination growth is likely to occur un- 
der mixed-mode loading. Therefore, a general formulation for decohesion elements 
dealing with mixed-mode delamination onset and propagation is also required. 
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Softening onset prediction Under pure Mode I, IT or III loading, the onset 
of damage at the interface can be determined simply by comparing the traction 
components with their respective allowables. However, under mixed-mode loading 
damage onset and the corresponding softening behavior may occur before any 
of the traction components involved reach their respective allowables, which is an 
issue that is usually neglected in the formulation of decohesion elements. Cui et al 
[37] have highlighted the importance of the interactions between interlaminar stress 
components when predicting delamination. It was shown that poor results are 
obtained using the maximum stress criterion. Therefore, a mixed-mode criterion 
accounting for the effect of the interaction of the traction components in the onset 
of delamination is proposed here. 

It is assumed that the initiation of the softening process can be predicted using 
the quadratic failure criterion [37], considering that compressive normal tractions 
do not affect delamination onset and using the operator defined in (19): 



This criterion has been successfully used to predict the onset of delamination 
in previous investigations [1], [2], [37]. 

The total mixed-mode relative displacement <S m is defined as: 

= \f#L + % + (<*> 3) 2 = \J tfshear + (^) 2 (28) 

where 8 s hear represents the norm of the vector defining the tangential relative 
displacements of the element. 

Using the same penalty stiffness in Modes I, II and III, the tractions before 
softening onset are: 

Ti = K8 h i = 1,2, 3 (29) 

Assuming S = T, the single-mode relative displacements at softening onset 
are: 


N 

«S = Ti (30) 
= = | (31) 

For an opening displacement £3 greater than zero, the mode mixity ratio 3 is 
defined as: 
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p = 8 j j^ (32) 

The mixed-mode relative displacement corresponding to the onset of softening, 
<5^, is obtained by substituting eqs. (28)-(32) into (27) and solving for <5 m , which 
gives: 


CO CO 

Mu 


St = < 


' 1+/? 2 

m + m ) 2 

S°shear <= $3 < 0 


6 * > 0 


(33) 


Clearly, pure mode loading is a particular case of the proposed formulation, as 
8° m = 8% for (3 = 0 (Mode I), and 8°^ = 8° shear for £3 = 0 (or when {3 — > 00 , Shear 
Mode) . 


Delamination propagation prediction The criteria used to predict delami- 
nation propagation under mixed-mode loading conditions are usually established 
in terms of the energy release rates and fracture toughness. There are established 
test methods to obtain the Mode I and II interlaminar fracture toughness. The 
Double Cantilever Beam Specimen (DCB) is used for Mode I. The End Notched 
Flexure (ENF) or the End Loaded Split (ELS) specimens are used for Mode II. For 
mixed-mode I and II, the Mixed-Mode Bending (MMB) test specimen is normally 
used. However, further research is required to assess the Mode III interlaminar 
fracture toughness, Gmc . Although some test methods have been suggested for 
the measurement of Mode III interlaminar fracture toughness, such as the Edge 
Crack Torsion [39] (ECT), there are important issues that need clarification, such 
as the determination of the transverse shear modulus G 23 , which is a parameter 
required for the analysis [14]. Furthermore, there is no reliable mixed- mode delam- 
ination failure criterion incorporating Mode III because there is no mixed-mode 
test method available incorporating Mode III loading. Therefore, most of the fail- 
ure criteria proposed for delamination growth were established for mixed-mode I 
and II loading only. For these reasons, and following Li’s work [4], [5], the concept 
of energy release rate related with shear loading, G shear = Gu + Gm, is used here. 

For mixed-mode loading the dependence of the fracture toughness on mode 
ratio must be accounted for in the formulation of decohesion elements. The relation 
between the mixed-mode interlaminar fracture toughness and the fracture surfaces 
of unidirectional laminates has been thoroughly examined using scanning electron 
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microscope analyses [21], [40]: for epoxy composites under pure Mode T loading the 
fracture surface is flat indicating cleavage fractures, whereas under pure Mode II 
loading the fracture surface exhibit hackles having an orientation of approximately 
45° with respect to the fiber direction. Under mixed- mode I and II the mechanisms 
are more complex, including both cleavage paths and hackles [21], [40]. 

The most widely used criterion to predict delamination propagation under 
mixed-mode loading, the power law criterion , is established in terms of an in- 
teraction between the energy release rates [41]: 


Gi 

Gic 


a 


+ 



(34) 


Reeder [21] performed mixed-mode bending (MMB) tests to measure the mixed- 
mode I and II interlaminar fracture toughness of composites, and obtained valuable 
experimental data to assess the several criteria proposed to predict delamination 
growth. The power law criterion obtained from (34) with a = 1 was found to be 
suited to predict failure of thermoplastic PEEK matrix composites because the 
results were comparable to the more sophisticated criteria, while using fewer inde- 
pendent variables. However, the power law criterion failed to accurately capture 
the dependence of the mixed-mode fracture toughness on the mode ratio occurring 
in epoxy composites using both a = 1 and a = 2. 

In order to accurately account for the variation of fracture toughness as a 
function of mode ratio in epoxy composites, the mixed-mode criterion proposed 
by Benzeggagh and Kenane [40] is used here ( B-K criterion). This criterion is 
expressed as a function of the Mode I and Mode II fracture toughness and a 
parameter p obtained from MMB tests at different mode ratios: 


Gic + (Guc — Gic ) ^ — Gc-, with Gt — Gi + Gu (35) 

If Mode III loading occurs the criterion is: 


Gic + (Guc — Gic ) ~q — ^ — Gc , with Gt — Gi + G shear (36) 

Figure 3 shows the predictions of the power law and B-K criteria for composites 
using a tough epoxy resin (IM7/977-2), a brittle epoxy resin (AS4/3501-6), and 
a thermoplastic resin (AS4/PEEK). The figure also includes the average of the 
experimental results obtained at different mode ratios (discrete points shown in 
Figure 3) and the a and p values used for each material. 
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G c (kJ/rrf) 



G„/G t 


Figure 3: Mixed- mode fracture toughness 


By using three parameters, Gic, Gnc and 77 , the B-K criterion represents 
the mixed-mode fracture toughness over a comprehensive range of mode mixity, 
whereas the two-parameter power law criterion that is usually implemented in 
decohesion elements can lead to inaccurate results over a large range of mode 
ratios. Furthermore, using the power law criteria with a = 2 in epoxy based 
composites, values of mixed-mode fracture toughness higher than Gnc occur in 
the interval 0.9 < Gh/Gt < 1.0. Experimental evidence (Figure 3) shows that 
the maximum mixed- mode fracture toughness occurs for Gjj/Gt = 1 , i.e., G JS ax = 
Gnc- Therefore, using of the power law criterion with a = 2 in epoxy based 
composites can lead to unconservative predictions in a small range of mode mixity 
ratios. 

Figure 3 shows that both the power law criterion using a = 1, and the B-K 
criterion provide reasonable results for the prediction of the mixed-mode fracture 
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toughness of AS4/PEEK composites. The power law criterion with a = 2 is 
clearly inadequate to predict the mixed- mode fracture toughness of AS4/PEEK 
composites. 

Based on the above results, the use of the B-K criterion in epoxy and thermo- 
plastic based composites is recommended. Taking into account that the maximum 
difference obtained by the application of the power law criterion using a = 1 to 
thermoplastic composites is 11.9% (at Gjj/Gt = 0.8), it is considered that the 
power law criterion using a = 1 can also be used, with the advantage of having 
one less variable than the B-K criterion. Therefore, both the B-K and the power 
law criterion are implemented in the decohesion element. 

The energy release rates corresponding to total decohesion are obtained from: 


G j 
Gu 
Gin 


L 

l 


6 3f 

v m 


T 'i,d 83 

(37) 

sU 


T2<18-2 

(38) 

Sm 


t idSi 

(39) 


Jo 


Using (24), (28) and (32) in equations (37)- (39) and substituting in (36) or in 
(34) the criterion for total decohesion can be established in terms of 8 m and (3 . 
Solving the equation for <5 m , the mixed-mode displacements corresponding to total 
decohesion, <5 j n , are obtained for the B-K criterion as: 


si = 


K8° 


Gic + (G no — Gic) 


0 


.2 \ V 


1 + /3 2 


83 > 0 


]/( 6{) 2 + ( 4 ) 2 ^ 83<0 

and for the power law criterion as: 

f 2(1 + /? 2 ) 


Si = 


K8Z 


1 \“ / y- 1 /** 


+ 


Gic ) V Guo J \ 

^/(< 5 f ) 2 + ( 4) 2 ^ 83 < 0 


83 > 0 


(40) 


(41) 
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Regardless of the criterion used, pure mode loading is a particular case of the 
proposed formulation, as 6^ = <5 t { for /3 = 0 (Mode T) and 8 ^ = 8{ hear for 6% = 0 
(or when (3 — > oo, Shear Mode). 

Constitutive equation for mixed-mode loading The constitutive equation 
for mixed-mode loading is defined by the penalty parameter K , the damage evolu- 
tion function d, the mixed-mode relative displacements corresponding to damage 
initiation and total decohesion, 8°^ and 8^ respectively, as: 


r s = D sr 8 r , with: 


D, r = < 


8 sr K <= 8%* < C 


(i - d)K + km. 3^1 <= r m < qr < € 


SaSar^S lK 


> sL 


d = 


6'JKt 1 - K.) 
- O 


, d € [0, 1] 


(42) 


(43) 


(44) 


It is worth noticing that Equation (43) avoids the interpenetration of the crack 
faces of the decohesion element for softening and fully open conditions. 

In order to define the loading and unloading conditions the state variable max- 
imum mixed-mode relative displacement , <5™ ax , and the loading function , F , are 
defined as: 


cm ax 

°m 

F(8m - C“) 


(trn ~ SZT) 
Sm ~ C* 


(45) 

(46) 


The mixed mode softening law presented above is a single- variable response 
similar to the bilinear single- mode law illustrated in Figure 2, defined by a damage 
evolution law (44), by the maximum mixed-mode relative displacement , (45), and 
by the loading function (46) . Only one state variable, the maximum relative dis- 
placement variable <5™ ax , is used to track the damage at the interface. By recording 
the highest value attained by <5 m , the unloading response is such as shown in Figure 
2. The relative displacements for initiation and ultimate failure are functions of 
the mode mixity /?, the material properties, and the penalty stiffness. 
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The mixed-mode softening law can be illustrated in a single three-dimensional 
map by representing Mode I on the Y-Z plane, and Shear Mode in the X-Z plane, 
as shown in Figure 4. 



Figure 4: Mixed-mode softening law 

The triangles 0 — N — 6$ and 0 — S — S'l hear are the bilinear response in Mode 
T and in shear mode respectively. In this three-dimensional map, any point on the 
0-X-Y plane represents a mixed- mode relative displacement. 

The map of all softening responses under mixed mode is illustrated in Figure 
5. The curve FI represents the tractions resulting from the displacements at the 
onset of damage given by (33), while the curve labeled G represents the ultimate 
relative displacements calculated with either (40) or (41). The triangle 0 — A — B 
is the bilinear softening law for a mixed-mode relative displacement of <5 m and the 
triangle 0 — C — D in Figure 5 is the Mode I bilinear softening response. It can also 
be observed that the effect of compression on the material response is neglected. 

2.3 Solution method for the non-linear problem 

The softening nature of the decohesion element constitutive equation causes diffi- 
culties in obtaining a converged solution for the non-linear problem. Furthermore, 
high penalty values can lead to large unbalanced forces and shoot the iteration 
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Tractions 



Figure 5: Map of softening response for mixed- mode delamination 

process beyond its radius of convergence. Crisfield et al [42] found that when 
using the Newton-Raphson method under load (with the arc-length method) or 
displacement control, the iterative solutions often failed to converge. In order to 
obtain convergence, a ’line search’ procedure with a negative step length was pro- 
posed. Other methods such as modified cylindrical arc-length method [43] and the 
nodal crack opening displacement (COD) control have been proposed [22], [30]. 

The Newton-Raphson method under displacement control and a ’line-search’ 
algorithm is used to solve the non-linear problem. Therefore, the decohesion ele- 
ment consistent tangent stiffness operator , must be calculated as: 

j^T dfki 

kizv 

U LL ZV 

As shown in Appendix A, the scalar components of the decohesion element 


(47) 
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consistent tangent stiffness matrix are: 


K L = B ri kKB rvz dr <= <5™*“ < 8° m (48) 

K L = K L,» + K L, <= C < < C where: (49) 

KT. = [ B sik S sr (1 - d)K + KSsz^- B rvz dT (50) 

Jr L — <->3 

K l,„ - ~ f B lit SrF(6 m - CTt 1 . 

Jr °m 

(51) 

(°m )(®m-U 

K L = J B Si Js^r^-KB rvz dT <= 5™ ax > Ha (52) 

where i,v are global degrees of freedom, fc, z represent node numbers and r,tc are 
related with the scalar components of the relative displacements. 

The function ^^,( 63 ) is defined as: 

W^) = Ssr (l - ^ 3 -hM) (1 _ ^ + M^) (53) 

The sum in (49) is a correction of the secant stiffness, K \ kivz , due to the 
damage growth occurring when F(6 m — <5™ ax ) = 1* Under unloading conditions, 
F(S m — <5^ ax ) = 0, and K kivz = K^ ki . It is worth noticing that the material 
tangent stiffness matrix included in K^ vz under single-mode unloading conditions 
corresponds to the secant stiffness shown in Figure 2, (1 — d)K. 

3 Simulation of delamination in fiber-reinforced 
composites 

The decohesion element proposed here was implemented in the ABAQUS Finite 
Element code [44] as a user- written element subroutine (UEL). To verify the ele- 
ment under different loading conditions, the double cantilever beam (DCB) test, 
the end notched flexure (ENF) test, and mixed-mode bending (MMB) tests are 
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simulated. The numerical predictions are compared with experimental data. The 
DCB test consists of pure Mode I delamination. The ENF tests measure pure 
Mode II interlaminar fracture toughness, and the MMB delaminate under Mixed 
Mode I and II. In the absence of Mode III loading, G s h ea r = Gn • To investigate 
the accuracy of the formulation in the simulation of delamination in different ma- 
terials, the DCB test is simulated for T300/977-2, a thermoset composite material, 
while the DCB, ENF and MMB simulations are conducted for PEEK/APC2, a 
thermoplastic matrix composite material. 

3.1 Mode I delamination growth for an epoxy composite 

The ASTM standard specimen used to determine the interlaminar fracture tough- 
ness in Mode I ( Gic ) is the double-cantilever beam (DCB) specimen [45] repre- 
sented in Figure 2 b). 

A DCB test specimen of a (0°)24, T300/977-2 carbon fiber-reinforced epoxy 
laminate, containing a thin insert at the mid-plane near the loaded end, is simu- 
lated. This specimen is 150-mm-long, 20 mm- wide, with two 1.98-mm-thick plies, 
and with an initial crack length of 55 mm. The DCB tests on this specimen were 
performed by Morais et al. and reported in [46]. The material properties are 
shown in Table 1. 


Table 1: Properties for T300/977-2 CFRP 


E ii 

#22 = #33 

CO 

0 

II 

CM 

O 

G 23 

150.0 GPa 

11.0 GPa 

6.0 GPa 

3.7 GPa 


^12 = ^13 

^23 

Gic 

T 

0.25 

0.45 

0.268 kJ/rri 2 

AbMPa 


In order to define the element constitutive equation, the penalty parameter 
and the interlaminar tensile and shear strengths must be determined. The choice 
of the penalty parameter can have an effect on the solution. Too low of a value 
leads to an inaccurate representation of the mechanical behavior of the interface, 
whereas high values can promote numerical errors related to computer precision. 
The optimum value for the penalty parameter is the largest value that does not 
lead to numerical problems. Some methodologies have been proposed to define the 
most adequate value for penalty parameters [47]. Based on previous investigations 
[36], a penalty K = 10 6 N/mm is used here. 
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Twenty-one node solid elements are used to simulate the DCB arms, and 18- 
node decohesion elements with a length of 1 mm are used along the interface. 

Figure 6 shows the results of 3 sets of experimental data and the numerical 
predictions. 

Load (N) 


60 

50 

40 

30 


20 


10 
0 

0 2 4 6 8 10 12 14 16 18 20 

Displacement (mm) 


+KV. 

V-v4-. 


T300/977-2 


fit 

/ 

ft-/ 


8 

a 

// 














Experimental 

Numerical 






Figure 6: Experimental and numerical results for T300/977-2 CFRP 

It can be seen that an excellent agreement between the experimental data and 
the numerical predictions is obtained. The averaged maximum load obtained in 
the experiments is 62.52 iV, whereas the maximum load predicted is 63.11 N . 
The unloading response is well reproduced by the numerical model, validating the 
unloading behavior of the constitutive equation proposed. 
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3.2 Mode I, Mode II and Mixed-Mode I and II delamina- 
tion growth for a PEEK composite 

The most widely used specimen for mixed-mode fracture is the mixed-mode bend- 
ing (MMB) specimen shown in Figure 7, which was proposed by Reeder and Crews 
[48] , [49] and later re-designed to minimize geometric nonlinearities [50] . This test 
method was recently standardized by the American Society for Testing and Mate- 
rials [51]. 




Loading arm 

-J— 


Specimen . 


Base 


V 




Figure 7: MMB test specimen 

The main advantages of the MMB test method are the possibility of using 
virtually the same specimen configuration as for Mode I tests, and the capability 
of obtaining different mixed mode ratios, ranging from pure Mode I to pure Mode 
II, by changing the length c of the loading lever shown in Figure 7. 

The 8-node decohesion element developed is used to simulate DCB, ENF and 
MMB tests in unidirectional AS4/PEEK carbon-fiber reinforced composite. The 
specimens simulated are 10 2- mm- long, 25.4 mm-wide, with two 1.56-mm-thick 
arms. The material properties are shown in Table 2, and a penalty stiffness K = 
10 6 N/mm is used. 
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Table 2: Properties for PEEK/APC2 


E n 

E 22 = E<x\ 

G 12 = Gn 

G 23 

V 12 = ^13 

122.7 GPa 

10.1 GPa 

5.5 GPa 

3.7 GPa 

0.25 


V23 

Gia 

Gnc 

T 

S 

0.45 

0.969 kJ/m 2 

1.719 kJ/rri 2 

80 MPa 

100 MPa 


The experimental tests were performed at different Gu/Gt ratios, ranging from 
pure Mode T loading to pure Mode II loading. The initial delamination length of the 
specimens (ao) and the Mixed-Mode fracture toughness obtained experimentally 
are shown in Table 3. 


Table 3: Experimental data 


Gu/Gt 

0% (DCB) 

20% 

50% 

80% 

100% (ENF) 

G c ( kJ/m 2 ) 

0.969 

1.103 

1.131 

1.376 

1.719 

a 0 (mm) 

32.9 

33.7 

34.1 

31.4 

39.3 


Models using 150 decohesion elements along the length of the specimens, and 4 
decohesion elements along the width, are created to simulate the ENF and MMB 
test cases. The initial size of the delamination is simulated by placing open de- 
cohesion elements along the length corresponding to the initial delamination of 
each specimen (see Table 3). These elements are capable of dealing with the con- 
tact conditions occurring for Mode II or Mixed-Mode I and II loading, therefore 
avoiding interpenetration of the delamination faces. The model of the DCB test 
specimen uses 102 decohesion elements along the length of the specimen. 

The different Gu/Gt ratios are simulated by applying different loads at the 
middle and at the end of the test specimen. The determination of the middle 
and end loads for each mode ratio is presented in Appendix B. The experimental 
results relate the load to the displacement of the point of application of the load P 
in the lever ( load-point displacement , Figure 7). Since the lever is not simulated, 
it is necessary to determine the load-point displacement from the displacement 
at the end and at the middle of the specimen, using the procedure described in 
Appendix C. 

The B-K parameter 77 = 2.284 is calculated by applying the least-squares fit 
procedure proposed in Appendix D to the experimental data shown in Table 3. 

Figure 8 shows the numerical predictions and the experimental data for all the 
test cases simulated, and Table 4 shows the comparison between the predicted and 
experimentally determined maximum loads. 
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Table 4: Experimental and numerical maximum loads 

Gh/Gt Pmax (experimental, N) P max (predicted, N) error (%) 
0% (DCB) 147.11 153.27 -4.2 

)% 108.09 86.95 19.6 

275.35 236.60 14.1 

518.66 479.86 7.5 

733.96 695.94 5.2 


100% (ENF) 





It can be concluded that a good agreement between the numerical predictions 
and the experimental results is obtained. The largest difference (19.6%) corre- 
sponds to the case of an MMB test specimen with Gu/Gt = 20%. This fact 
is not surprising, since the largest difference between the fracture toughness ex- 
perimentally measured and the one predicted using the B-K criterion occurs for 
Gii/Gt = 20% (see Appendix D). 

4 Conclusions 

A method for the simulation of progressive delamination based on decohesion el- 
ements was presented. Decohesion elements are placed between layers of solid 
elements that open in response to the loading situation. The onset of damage and 
the growth of delamination can be simulated without previous knowledge about 
the location, the size, or the direction of propagation of the delaminations. A 
softening law for mixed-mode delamination that can be applied to any interaction 
criterion was proposed. The criterion uses a single state variable, the maximum 
relative displacement, to track the damage at the interface under general loading 
conditions. For the linear and quadratic power law criteria, the material properties 
required to define the element constitutive equation are the interlaminar fracture 
toughnesses and the corresponding strengths. The B-K interaction law requires 
additionally a material parameter r] that is determined from standard mixed-mode 
tests. 

Three examples were presented that test the accuracy of the method. Simula- 
tions of the DCB and ENF test represent cases of single-mode delamination. The 
MMB tests that were simulated have various proportions of Mode I and Mode II 
loading conditions, and were simulated using a three-parameter criterion for de- 
lamination propagation. The examples analyzed are in good agreement with the 
test results and they indicate that the proposed mixed-mode criteria can predict 
the strength of composite structures that exhibit progressive delamination. 
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Appendix A: Determination of the consistent 
tangent stiffness matrix 

The tangent stiffness matrix used in solution of the non-linear problem is de- 
fined by: 


_ dfki 


du zv 

and the decohesion element internal is load vector is: 


(54) 


fki J B S ifcT s dV J B si kD sr B rvz dT u z 


(55) 


The consistent tangent stiffness matrix is obtained for the different states of 
the interface (no damage, softening or open interface): 


Case 1 No damage, <5™ ax < 6 ^ and D sr = 8 sr K. 


From (54) and (55): 


K Lv = J B rikKB rvz dV 


(56) 


Case 2 Softening, 8° m < 5“ ax < 6^ and D sr = 6 sr 
From (54) and (55): 


(1 -QK + KdS*^- 


= J B sik D sr B rvz dT + J B si k B rvz dT u zv 


K 


T 


f dD sr „ 

B S ikD sr B rvz dT + / B S ik— 8 r dT 


du zq 


K t = KT + KT 

kizv L kizv kizv 


The definition of the term in Kl is obtained as follows: 

OUzv ^ kizv 


(57) 

(58) 

(59) 


dDsr 

du zv 

dd 

du„, 


dD sr dd 
dd du zv 
dd dS w 

dSyj dll Z y 


= -KSsr 1 - <5 


sr \ x ^s3‘ c 
—0 3 

J^dVT' , R 
d6T x dS„ 1 wzv 


(Ss)\ dd 




= -A'KfhS 


dd 

du zv 


(60) 

(61) 
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From (44): 


3i SUC 

<?«"“ - to 


From (28): 


(62) 


dS w 


F($m - - 6*. + ^-8, w ) = F(8 m - qr&VZQs) 

Um ^3 


(63) 


Using (62) and (63) in (60): 


dDsr 

dtizv 


-K**(8 a )— 
3 (6" 


sLst 


‘WL - K 


■F(s m - «;r) 5 r'C(fe)B. 

Om 


(64) 


From (64) and (58)- (59): 


*L. 


- [ B sik 8 r F(8 m -8™ x )^. 

Jr 

K6LSL 


Kk 


\8Z*)\8l-8 

- f B sik 8 r F(8 m - C“) 


9l(6 3 )*Z(6 3 )B wgv dr 


KSLSt 


" 4 ^ srw (83) B wzv dF 


'(S£*nsL-K 

Case 3 Open interface, <§™ ax > 8^ m and D sr = 8 s ^8^ r ^^- K. 
From (54) and (55): 


K 1 


B s ik 8 s 3 8% r ~ — T^~ K B rvz dr 
' — Os 


(65) 


(66) 


(67) 


32 



Appendix B: Determination of middle and end 
load in MMB tests 

The length of the lever used in the MMB test, c, was obtained taking into 
account the weight of the lever. Since the lever is not simulated in the numer- 
ical models, the lengths corresponding to the different mode ratios need to be 
calculated. 

The mode mixity ratio, ft, is defined as a function of the energy release rates 
as: 


G 


ft = 


ii 


Gi 

Gn 


(i-«) 


( 68 ) 


Gj + Gjj Cr a ft 

The relation between Mode I and Mode IT energy release rates in a MMB test 
specimen is [52], [49]: 


Gi 

Gu 


3c- zy f ^ 1 
7TT) ' forc£ 3 


(69) 


From (68) and (69), the length of the lever can be obtained as a function of 
the mode mixity ratio and specimen length (Z): 


c = 


(H/3(¥)+l 


(70) 


Table 5 shows the lengths of the lever for each mode mixity use in the numerical 
model and in the experiments. 


Table 5: Lengths of the lever 


Gnj Gt 

20% 

50% 

80% 

c (numerical, mm) 

109.4 

44.4 

28.4 

c (experimental, mm) 

97.4 

42.2 

27.6 


Neglecting the weight of the lever, the middle and end loads, P m and P e re- 
spectively, are obtained as a function of the total load P as [49] : 


P n = P(^) (71) 

P. = Pj (72) 
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Therefore: 


From (73) and (70): 


Pm _ c + l 

Pe C 


( 73 ) 


Pm _ o 6k + a/3k(1 - k) 

Pe 3 + 9/C + 8y / 3/«(l — k) 


(74) 


Using (74) the relation between the load at the middle and at the end of the 
specimen is obtained for the different mode ratios: 


Table 6: Ratio between middle and end loads 


Guj Gt 

20% 

50% 

80% 

Pm/ Pe 

1.46 

2.14 

2.79 


The pure mode load components have been shown [49] to be given as a function 
of the total load applied to the specimen by: 


Pi = P 
Pn = P 



From (70), (75) and (76): 


(75) 

(76) 


(l - K + 2 ^3 K (1 - «)) P 
13 K — 1 

g (6 k + \/3k (1 — k)^ P 

P T 1 — 

3 13 ft — 1 

The Mode T and Mode II loads for each case are shown in Table 7. 


(77) 

(78) 


Table 7: Mode I and Mode II load components 


Gu/ G t 

20% 

50% 

80% 

Pi 

1.37P 

0.41P 

0.17P 

Pii 

3.15P 

1.87P 

1.56P 
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Appendix C: Determination of the load-point 
displacement 

The information available from the MMB test relates the load to the displace- 
ment of the load-point (Figure 7). Since the lever is not simulated in the numerical 
models, it is necessary to calculate the load-point displacement using the informa- 
tion available from the numerical models of the MMB test specimens. 

The load-point displacement, A^p, is obtained from the pure mode displace- 
ment components, A / and A//, as [53]: 

Alp = A r) A; + (A) A,J (79) 

Using simple beam theory, Mi and Crisfield [52] calculated the Mode II dis- 
placement component, A//, as a function of the displacement at the middle of the 
MMB test specimen, Am- 


Therefore: 


A n — Am + -A/, for a < l 


(80) 


Alp = ^ ^ A/ + ^ ^ ^Am + -A/^ , for a < l (81) 

The values of A/ and Am are computed from the numerical model. The load- 
point displacement is then calculated using (81), being this procedure valid only 
for crack lengths (a) smaller than half the length of the specimen (/). 
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Appendix D: Experimental determination of the 
parameter v for the B-K delamination criterion 


The problem consists in determining 77 from a set of experimental data using 
the polynomial: 


p(-^A) — Gic + (Guc - Gfc) ( -^A) (82) 

A least-square fit is proposed. Considering the pair ( (^V ) » ( Gt)jJ as the 

experimental data and n as the number of data points, the problem can be posed 
as the minimization of: 


S=E 

3 = 1 

Considering ^ = 0 : 



Gic — (Guc — Gjc ) 



j 


(83) 


n 


E 

3 = 1 


( G t)j 


Gic — {Guc — Gic ) 



v 

j 



0 (84) 


The value of 77 is then obtained from the solution of equation (84). For the 
experimental data used (Table 3), 77 = 2.284. The application of the B-K criterion 
over the entire range of mode-mixity ratios is shown in Figure 9. 

It should be noticed that only one experimental point was available for each 
loading condition, and that the largest difference between the B-K criterion and 
the experimental results occurs for Gu/Gt = 0 . 2 . 
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